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1. Introduction 



The problem we are addressing is the development of thermodynamic ocean 
models which will simulate the horizontal variability of vertical mixing in 
the upper ocean. These models will be capable of simulating the upper ocean 
response to local atmospheric forcing on synoptic to seasonal time scales. 

The spatial response of the models will be consistent with resolution of the 
atmospheric fields used as forcing; thus the atmospheric fields should accur- 
ately represent the mesoscale to large scale variability of the atmosphere. 

The fields which are needed to force the models include the solar radiation, 
the total heat flux and the wind speed at the ocean surface. The atmospheric 
prediction models and the analyses of atmospheric observations used by Fleet 
Numerical Weather Central (FNWC) to produce operational weather forecasts can 
provide the atmospheric forcing for these ocean models. 

This report describes a software system which has been designed to produce 
the time series of atmospheric fields needed to drive a one-dimensional, tur- 
bulence closure, thermodynamic ocean model (Garwood, 1977a), from the FNWC 
analysis of the surface marine wind field, and the FNWC predictions of solar 
radiation, and total heat flux at the ocean surface (FNWC, 1974). Steve Pazan, 
the North Pacific Experiment (NORPAX) data manager, has provided us with the 
FNWC fields for the Anomaly Dynamics Study (ADS) region (30-50N, 140E-130W) , 

Fig. 1. This data set is on a 11 x 41 grid for the years 1976 and 1977 (Anomaly 
Dynamics Study, 1978). Since March 1978, we have been receiving the fields 
directly from FNWC on the 63 x 63 polar stereographic Northern Hemispheric grid. 
Fig. 2. The results discussed in this report have been obtained using the ADS 
data set; however, the software system is directly applicable to the FNWC data 
set. 

The diurnal heating of the upper ocean and its coupling to the synoptic 
scale wind field have been shown to be important in predicting the seasonal 



and annual cycles in the upper ocean (Garwood, 1977b) due to the nonlinearity 
of the ocean response to the atmospheric forcing. The mixed layer model must 
be supplied with atmospheric forcing at least every 3 hours to properly re- 
solve the ocean response to the diurnal signal. The sea level (10 m) winds, 
from the FNWC analyses, are saved at 6 hour intervals. Instantaneous values 
of the surface heat fluxes are recorded every 12 hours (every 6 hours since 
March 1978). These fields are permanently archived on magnetic tape and are 
available to the scientific community. The mission of FNWC is to provide 
environmental support for fleet units; hence, the archiving of fields for later 
research has not been a primary consideration. As a result all the FNWC fields 
were not archived for every 6 or 12 hour interval. In addition, some archived 
fields were duplicated when test runs of the models or analyses were inadver- 
tently recorded. The number of missing or duplicate fields is normally less 
than 10% of the total for any month. The incidence of missing or duplicate 
fields is becoming less frequent as the interest in using the archived fields 
increases. Nevertheless, part of the function of our data management system 
must be to locate and correct missing and duplicate fields. Additionally, it 
must be noted that while the surface fluxes are of primary importance to ocean 
modeling, they are of secondary importance in atmospheric models. Thus the 
FNWC fields must be carefully screened for any bias or systematic error which, 
while not important to the atmospheric model in operational integrations, could 
overwhelm the ocean model during longer term integrations. 

2. System Description 

The software system, which converts the FNWC atmospheric boundary layer 
fields into the forcing required for the Garwood mixed layer model, consists 
of nine programs divided into three quasi-independent subsystems, Fig. 3. 

This approach provides natural breaks in the data management and gives the 
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system greater flexibility by limiting the amount of computer time required 
for any one run. This method also allows the system to be used with alter- 
native atmospheric data sets and to accomplish a variety of functions which 
can easily be altered to accommodate the different data sets. 

The first subsystem, REFORMAT, reformats, edits and sorts the FNWC fields, 
Fig. 3a. The resulting data set is compatible with the IBM 360/67 computer 
system at the Naval Postgraduate School, and can be read more than an order of 
magnitude faster than the original. The next subsystem, RETRIEVE, retrieves 
any desired subset in space and time of the data set generated by REFORMAT, 
and interpolates values for the missing fields which were flagged by REFORMAT, 
Fig. 3b. The final subsystem, INTERPOLATE, interpolates the FNWC fields to 
the time interval required for the model. Fig. 3a. Some of the fields require 
specialized interpolation schemes due to the distinctive character of the 
fields. 

Each subsystem will be discussed in detail in the following sections. All 
the programs in this system were written in FORTRAN IV with the exception of 
the sorting program, which is an IBM utility program. All the data files are 
written in binary format, which yields a considerable savings in file space 
and processing time. Each data file contains a header record, which specifies 
the space and time domain of the data, followed by the data records. 

2a. Reformat Subsystem 

The REFORMAT subsystem consists of three programs that transform the FNWC 
fields into the format used by our software system. The first of these programs 
RFMT, unpacks the FNWC or ADS data tapes, and converts the data fields to 4-byte 
real variables with the proper physical units. The FNWC fields are archived 
as 4 character octal variables to conserve space on the data tapes. The 
second program, EDIT, tests the RFMT generated data file for missing and 
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duplicate records. Each data record contains an entire grid of data points 
for one surface parameter at a given time. Duplicate records are dropped 
from the file and missing records flagged by inserting a grid of 999.9 values 
for the missing record. The data file generated by the EDIT program contains 
a record for each surface parameter for every time interval. The 999.9 grids 
which have been inserted for missing records can be replaced with reasonable 
values or deleted from the file later in the processing. This arrangement 
allows considerable flexibility in handling the data. The final program in 
the subsystem, SORT, arranges the data records in chronological order. The 
data files, which results from the REFORMAT subsystem, contains a data record 
for every FNWC field at each time interval. This data file can be read at 
least 20 times faster than the original FNWC or ADS tapes. 

2b. Retrieve Subsystem 

The RETRIEVE subsystem consists of two programs which can extract any 
subset of the files generated by the REFORMAT subsystem and replace the missing 
(999.9) records with useable data values. The first program, RTRV, extracts the 
desired data as specified by a control card. The control card designates the 
initial and final latitudes and longitudes, the date and time of the first 
record and the number of records to be extracted. The requested subset is 
retrieved and a data file is created for each of the FNWC surface parameters. 

The second program, CRCT, searches the files generated by RTRV and replaces 
the missing (999.9) fields with data values derived by linear interpolation. 

The RETRIEVE subsystem generates four data files; one each for the zonal surfac 
wind component, the meridional surface wind component, the surface solar radi- 
ation and the total heat flux at the surface. Each data file contains a com- 
plete time series for a FNWC field. These can now be interpolated to the 
time step of the ocean model. 
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2c. Interpolation Subsystem 

The INTERPOLATION subsystem is composed of four programs which interpolate 
the FNWC fields to 1 hour intervals for the Garwood mixed layer model. The 
characteristics and archiving frequency of the individual FNWC fields must be 
considered in choosing an interpolation scheme for each field. For example, 
the wind fields can be interpolated using a straightforward cubic spline 
method; however, the interpolation of the solar radiation requires a more 
complicated procedure. The cubic spline routines are taken from the Inter- 
national Mathematics and Statistics Library (IMSL, 1979). The interpolation 
to hourly intervals is a two step process. In the first step the coefficients 
of the cubic spline are calculated at each knot, which is defined as an inter- 
polated point. In the second step, the coefficients are used to evaluate the 
cubic spline between the knots. 

The wind components are interpolated using the program FRCFI, which serially 
interpolates each grid point and outputs a time series of hourly records; 
each record containing all the grid points for that hour. After FRCFI is run 
for both wind components, the program, WNDSP, is used to calculate the wind 
speed at each grid point from the wind components. A typical wind speed 
profile is shown in Fig. 4. The wind speeds computed from the original 6 hour 
wind components are marked with a vertical bar. 

The program, A11I, is used to interpolate the solar radiation to 1 hour 
intervals. The solar radiation has a semi-diurnal periodicity and only one 
instantaneous value each 12 hours, at OOGMT and 12GMT; thus a simple inter- 
polation scheme is not feasible. Instead we make use of Malankovich* s formula 
to estimate hourly values of the solar radiation. Malankovich 1 s formula is 
given as 
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where 1^ = surface solar radiation of the i-th hour from local noon, 

I = surface solar radiation at local noon, 
n 

q = water vapor absorptance for the atmosphere which is taken as 0.7, 
= sin $ sin 5 + cos <f) cos 5 cosh^, . > 

(() = latitude, 

5 = solar declination, 

h i = hour angle for the i-th hour (15°/hour). 



Values of y— for 6 = 23° (spring equinox) ,5=0° (summer solstice) and 6 = -23° 
n 

(vernal equinox) at 40N are shown in Fig. 5. Malankovich’ s formula is used to 
calculate the solar radiation each hour from one hour after sunrise to one hour 
before sunset based on the 12 hour FNWC value which falls closest to local 
solar noon. This method of calculating the hourly solar radiation assumes that 
the atmospheric water vapor and cloudiness are constant during the hours of 
daylight. Thus high frequency variability in the instantaneous measurements 
of the solar flux due to clouds and water vapor are not included. Rather the 
assumption is that the significant time scales for heating the ocean are related 
to the synoptic storm systems and the diurnal cycle of heating and cooling. 
Accumulation of these synoptic and diurnal scale heat fluxes result in the 
synoptic to seasonal scale variability in the heat content of the upper ocean 
through vertical mixing. Fig. 4 is a typical example of the interpolated 
solar radiation, with the FNWC values shown as f, x". The program, A18I, inter- 
polates the 12 hour FNWC fields of the total heat flux to hourly intervals. 

This interpolation is a three step process because the solar flux must be 
handled separately. First, the original 12 hour solar radiation values are 
subtracted from the corresponding total heat fluxes to form a residual heat 
flux. The residual heat flux includes the sensible, latent and long wave heat 
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fluxes; each of which is a smoothly varying function with synoptic scale 
variability. Therefore we can interpolate hourly values of the residual heat 
flux from the 12 hour values. The method employed is identical to that used 
for the wind components. The last step in the procedure is to add the hourly 
solar radiation fields, created by All I, to the hourly residual heat fluxes 
to form hourly total heat flux fields. A typical result is shown in Fig. 4. 

3. Conclusions 

This system accomplishes the task of manipulating the FNWC atmospheric 
boundary layer fields into time series of forcing for the Garwood mixed layer 
model. The divisions of the system into the particular subsystems and programs 
discussed above facilitates data control during both system testing and pro- 
duction. It also allows the system to be easily restructured to accommodate 
different data sets, or different methods of data manipulation. The compart- 
mentalization of the system simplifies its use by other groups; few if any 
changes would have to be made and these would involve a minimum number of 
programs or subsystems. A brief documentation of the programs, and auxiliary 
programs and subroutines is given in the Appendixes. 

Preliminary tests (Elsberry, Gallacher and Garwood, 1979) indicate that 
the system accurately handles the FNWC fields and prepares useable forcing for 
the Garwood model. Testing is being continued to determine the feasibility 
of long term ocean prediction using the FNWC fields. We need to consider both 
the accuracy of the low frequency components and the lack of high frequency (>2 cpd) 
resolution in the FNWC fields. The surface heat flux fields are of primary 
importance to ocean modeling; however, they produce only second order effects 
in the atmospheric models. Thus a bias could exist in the surface heat fluxes that 
may not impair the atmospheric forecasts but which could result in systematic 
errors in long term ocean simulations. Additionally, we know that during the 
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spring transition and the fall deepening of the ocean mixed layer, the phase 



difference between the wind forcing and the surface heat fluxes can critically 
affect the ocean response on synoptic time scales (Elsberry and Garwood, 1978). 
It is not known if the high frequency variability is also important for sea- 
sonal and annual simulations. The answer to this has serious implications 
since the FNWC fields can not resolve diurnal and higher frequency variability. 
Also the peak values tend to be reduced by the objective analysis schemes 
which are used to produce the wind speed and surface heat flux fields. Several 

tests and experiments are being conducted to resolve these questions. 

As we have discussed, the software system can handle ADS or FNWC generated 
surface fields. The system can be adapted, with minimum modifications, to 
handle satellite data or surface data from ocean experiments such as MILE, 
FRONTS, and STREX. Thus we can look at the results from a variety of sources 
on several space and time scales. 
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Appendix A 



PROGRAM: 


RFMT 


PURPOSE: 


To reformat the FNWC or ADS data tapes, creating a new data set 
which is compatible with the IBM 360/67 computer system at the 
Naval Postgraduate School. 


INPUT: 


ADS or FNWC data tape. 


OUTPUT: 


A data file in binary format with all data values in physical units 
and unnecessary header information deleted. 


PROCESS: 


The header record and data records are read from the ADS or FNWCS 
tapes. The scale factor given in the header record is used to 
convert the data elements from 4 character octal constants to real 
*4 variables. A shortened version of the header record is output 
followed by the data records. All the output records are binary 
(i.e. written using FORTRAN unformatted writes statements). 





Appendix B 


PROGRAM: 


EDIT 


PURPOSE: 


To remove duplicate records and fill in missing records in the 
data file. 


INPUT: 


Data file generated by RFMT. 


OUTPUT : 


A binary data file with no duplicate fields and a data record for 
every time interval and a printed summary of the actions taken by 
the program. 


PROCESS: 


The input file is scanned for duplicate date-time fields and for 
missing date-time groups. A duplicate date- time group is defined, 
in the program, as a date-time group which is less than the date- 
time group which is expected for the data record which the program 
has read. When a duplicate record is located the catalog number 
and date-time group are printed on the summary sheet and the data 
record is not written to the output file. A missing date-time 
group is located when the date-time group of the data record is 
greater than the date-time group expected by the program. When 
a missing record is located, the missing date-time group and 
catalog number are printed on the summary sheet and a data record 
filled with 999.9 values is written to the output file. 
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Appendix C 



PROGRAM: 


SORT 


PURPOSE: 


To put the data file generated by EDIT in chronological order. 


INPUT: 


Data file generated by EDIT. 


OUTPUT : 


A data file in chronological order. 


PROCESS : 


The input file is sorted into ascending order on the date-time 
group by the IBM utility sort program. 





Appendix D 


PROGRAM: 


RTRV 


PURPOSE: 


To retrieve from the complete data set of forcing functions a time 
series for each forcing function at the grid points specified on the 
control card. The time series starts on the day specified and con- 
tains the number of days specified by the programmer. 


INPUT: 1. 


ADS forcing tape which contains the wind velocity components at 6 
hour intervals, the total heat flux, the solar radiation, the 
evaporative heat flux, the cumulative precipitation and the total 
cloud cover every 12 hours. 


2. 


Control card containing the initial and final latitude and longitude 
of the grid point (s) to be retrieved, the date-time of the beginning 
of the time series and the number of days to be extracted. 


OUTPUT : 1 . 


Time series of the solar radiation (All) for the grid point (s), 
starting time and number of days specified in the control card. 


2. 

3. 

4. 


Same as (1) but for the total heat flux (A18) 

Same as (1) but for the E-W wind component (A29) 

Same as (1) but for the N-S wind component (A30) 


PROCESS: 


The program reads the control card and translates the latitude and 
longitude into array indices. The beginning date-time is translated 
into Julian date format and an ending date-time is computed using the 
number of days to be retrieved. The input file is searched and all 
records with catalog numbers All, A18, A29, A30 which fall between 
the starting and stopping time are extracted and placed in the output 
files by catalog number. 
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Appendix E 



PROGRAM: 

PURPOSE: 

INPUT: 

OUTPUT: 

PROCESS: 



PROGRAM: 

PURPOSE: 

INPUT: 

OUTPUT: 

PROCESS: 



CRCT 

To replace any missing fields in the files created by RTRV. 

Any one of the files generated by RTRV. 

A new file which contains a useable value for every date-time. 

The input file is scanned for records containing a data value of 
999.0. When a record containing this value is found the 999.0 value 
is replaced with a value derived by linear interpolation using the 
two nearest values in time. The 999.0 values were inserted for all 
fields which were missing in the original data tapes by the program 
EDITA. The program outputs a new time series which consists of all 
the data records from the input time series with the 999.0 values 
replaced with interpolated values. Note: This program must be run 

for each file generated by RTRV1D. 



Appendix F 



All I 

To interpolate the values of solar radiation to 1 hour time intervals. 
Time series of solar radiation (All) generated by CRCT. 

Time series of solar radiation at 1 hour intervals. 

One day of data points are read from the input file, sunrise and sun- 
set are computed for that day. The data value which is closest to 
local noon is located. Data values are then computed for every hour 
between sunrise and sunset using Malankovich 1 s formula, which states 
that 

1 _ _1 

(sina, sina ) 

1(h) = sina^ I(n) q n 

sina 

n 

where sina^ = sin <f> sin 6 + cos $ cos 6 cos h 
<f> = latitude 
6 = solar declination 
h = hour angle = (TIME-NOON) *15° /hr 
q = 0.7 

I(n) = data value nearest noon. 

The interpolated values are then written to a file with the values 
for night time hours set to zero. 
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Appendix G 



PROGRAM: 


A18I 


PURPOSE: 


To interpolate the values of total heat flux to 1 hour intervals. 


INPUT: 1. 
2. 
3. 


Time series of total heat flux (A18) generated by CRCT. 
Time series of solar radiation (All) generated by CRCT. 
Time series of solar radiation (All) generated by A11I. 


OUTPUT: 


Time series of the total heat flux at 1 hour intervals. 


PROCESS: 


The 12 hour interval values of the solar radiation are subtracted 
from the 12 hour values of the total heat flux to obtain the residual 
heat flux. This residual heat flux is interpolated to 1 hour time 
intervals using cubic splines (the IMSL cubic spline routines are 
used). The values of solar radiation at 1 hour intervals are now 
added to the interpolated residual heat flux to obtain the total heat 
flux at 1 hour intervals. 



Appendix H 



PROGRAM: 


FRCFI 


PURPOSE: 


To interpolate any field which does not require special handling to 
a 1 hour time interval. 


INPUT: 


Time series of field to be interpolated (in our case this is one of 
the wind components). 


OUTPUT : 


Time series of forcing at 1 hour time intervals. 


PROCESS: 


The input time series is read and interpolated to 1 hour intervals 
using the IMSL cubic spline routines. The results are output to a 
new file. 



Appendix I 



PROGRAM: 


WNDS 


PURPOSE: 


To form a wind speed (A27) from the E-W and N-S wind components. 


INPUT: 1. 

2. 


The N-S wind components (A30) at 1 hour intervals. 
The E-W wind components (A29) at 1 hour intervals. 


OUTPUT : 


The wind speed at 1 hour time intervals. 


PROCESS: 


The wind speed is computed from the wind components (i.e. 

W = (W^ + Vr)^). The results are then output to a new file, 
x y J ' 
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Appendix J 



PROGRAM: 


TSERS 


PURPOSE : 


To plot time series of the forcing functions used by the model. 


INPUT: 1. 

2. 

3. 

4. 


The solar radiation (All) time series at 1 hour intervals. 

The total heat flux (A18) time series at 1 hour intervals. 

The wind speed (A27) time series at 1 hour intervals. 

A control card giving the min and max values of the x and y variables 

for each plot. 


OUTPUT: 


A versatec plot of the solar radiation, total heat flux and wind 
speed as functions of time. 


PROCESS: 


The program inputs the control card and time series and uses versatec 
plot calls to build a complete plot with axes, labels and titles. 




Appendix K 



SUBROUTINE: LOCADS 



PURPOSE: 


To convert from latitude (s) and longitude (s) to array indices. 


INPUT: 


ALATC1 - starting latitude of the grid point (s) of interest. 

AL0NC1 - starting longitude of the grid point (s) of interest. 

ALATC2 - ending latitude of the grid points of interest. 

AL0NC2 - ending longitude of the grid points of interest. 

DELAT - displacement in degrees latitude between desired grid points. 
DELON - displacement in degrees longitude between desired grid points 


OUTPUT: 


ILAT1 - array index of first grid point for latitude index. 

IL0N1 - array index of first grid point for longitude index. 

ILAT2 - array index of last grid point for latitude index. 

IL0N2 - array index of last grid point for longitude index. 

LATD - increment between desired grid points for latitude index. 
LOND - increment between desired grid points for longitude index. 


PROCESS: 


The subroutine computes the starting and ending values of the array 
indices and the increment for the indices based on the latitudes, 
longitudes and displacements passed to the subroutine. For the one 
dimensional case (only one grid point desired) ALATC1 and AL0NC1 are 
set to the latitude and longitude of the grid point. ALATC2 = AL0NC2 
= DELAT = DELON = 0 (or blank). In this case the array indices of 
the grid point are returned from the subroutine in ILAT1 and IL0N1. 

For the two dimensional case, ALATC1 and AL0NC1 are the latitude 
and longitude of the lower left hand comer (point furthest to the 
south and west) of the set of grid points desired; ALATC2 and AL0NC2 
are the latitude and longitude of the upper right hand corner of the 
set of grid points; DELAT is the latitude displacement between desired 
grid points and DELON is the longitude displacement between desired 
grid points. If DELAT = DELON = 0.0, the subroutine sets ILATD and 
ILOND to one so that no grid points are skipped. If DELAT and DELON 
were set to twice the spacing between the original grid points, ILATD 
and ILOND are set to 2 so that every other grid point will be used 
by the calling program. 
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Appendix L 



FUNCTION: JULIAN 



PURPOSE: 


To convert from a calendar date-time to a Julian date-time. 


INPUT: 


Calendar date-time in the form YYMMDDTT 
where YY = year - 1900 i.e. 78 = 1978 

MM = month 
DD = day 

TT = time in hours i.e. 1200 = 12 


OUTPUT: 


Julian date^time group in the form YYDDDTT 
where YY is the same as above 

DDD is the day of the year i.e. Jan 10 = 010 
TT is the same as above. 


PROCESS: 


The function calculates the Julian date-time group which corresponds 
to the calendar date-time group sent to the function. Leap years 
are taken into account. 

Appendix M 


FUNCTION: 


ICLNDR 


PURPOSE: 


To convert from a Julian group in the form YYDDDTT 
(see JULIAN for definitions). 


OUTPUT: 


Calendar date-time group in the form YYMMDDTT 
(see JULIAN for definitions). 


PROCESS: 


The function calculates the calendar date-time group which corresponds 
to the Julian date-time group sent to the function. 

Appendix N 



SUBROUTINE: SUN 



PURPOSE: 


To calculate sunrise in GMT and the number of hours of daylight 
for any latitude, longitude and day of the year. 


INPUT: 


ALAT - latitude of grid point 
ALON - longitude of grid point 
DAYNO - day of the year. 


OUTPUT: 


DEC - solar declination 

ALATR - latitude of point in radians 

SUNR - time of sunrise in GMT 

DLHRS - number of hours between sunrise and sunset 
IGDAY - 0 if local DAYNO and GM DAYNO of sunrise are the same 
-1 if GM DAYNO of sunrise is less than local DAYNO 
+1 if GM DAYNO of sunrise is greater than local DAYNO, 


PROCESS: 


The subroutine calculates sunrise and the number of hours of daylight 
for any latitude and longitude in the Northern Hemisphere for any day 
of the year. Sunrise and sunset will be symmetric with respect to 
local noon. 
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rig. 1 The Anomaly Dynamic Study (ADS) region (10-50N, 140r:-130U) . 




Fig. 2 The 63 x 63 polar stereographic grid for the Northern Hemisphere used 
• by FNWC. 
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REFORMAT 

SUBSYSTEM 




Fig. 3a The REFORMAT subsystem consists of three programs: RFMT, which refor- 

mats the data; EDIT, which checks for missing and duplicate fields; 
and SORT, which resequences the data set. 
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RETRIEVE 

SUBSYSTEM 




Fig. 3b The -RETRIEVE subsystem consists of two programs: RTRV, which retrieves 

a subset of the complete FNWC or ADS grid; and CRCT, which interpolates 
values to fill missing fields. 
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interpolate 

SUBSYSTEM 




Fig. 3c The INTERPOLATE subsystem consists of four programs: A11I, which 

interpolates the solar radiation; A18I, which interpolates the total 
heat flux; FRCFI, which interpolates the wind components; and WNDSP, 
which computes the wind speed from the vector components. 
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Fig. 4 The interpolated wind speed, total heat flux and solar radiation at the 
surface for 32N, 175W. 
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Fig. 5 Values of computed by Malankovich’ s formula for 40N, 6 = 23°, 0°, -23' 

n 
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